!******************************************************************************************
!Subroutine amoeba from Numerical Recipes (NR) for Fortran 90
!******************************************************************************************
!------------------------------------------------------------------------------------------ 
! NOTE: This subroutine performs the Nelder-Mead downhill Simplex method. This 
! multidimensional minimization method does not require derivatives and relies only 
! on function evaluations. This program code is based on NUMERICAL RECIPES (p. 404-405). 
!It is slightly modified to allow for non-uniform changes across different parameters.
! The method consists in minimizing the function FUNC in N dimensions by
! the downhill simplex method of Nelder and Mead. The (N+1)*N matrix P is input (N is the
! # of parameters to be estimated). Its N+1 rows are N-dimensional vectors which are the
! vertices of the starting simplex. Also input is the vector Y of length N+1, whose 
! components must be pre-initialized to the values of FUNC (here the likelihood function)
! evaluated at the N+1 vertices (rows) of P, i.e., the value of FUNC at the 
! vertices of the starting simplex. FTOL is the fractional convergence tolerance 
! to be achieved in the FUNCTION VALUE. On output, P and Y will have been reset to N+1 
! new points all within FTOL of a minimum function value, and ITER provides the number 
! of iterations (function evaluations) taken.
!------------------------------------------------------------------------------------------
SUBROUTINE amoeba(p,y,ftol,iter)


USE mod_param_initial
USE mod_nrutil
USE mod_parameters
USE mod_types


IMPLICIT NONE


INTEGER, INTENT(OUT):: iter
REAL(dp), INTENT(IN):: ftol 
REAL(dp), DIMENSION(num_param+1), INTENT(INOUT):: y
REAL(dp), DIMENSION(CAP_N):: logfct
REAL(dp), DIMENSION(num_param+1,num_param), INTENT(INOUT):: p
INTEGER, PARAMETER:: itmax = 20000000
!in NR it is set to 1.0d-10
REAL(dp), PARAMETER:: tiny = 1.0d-10 
INTEGER:: ihi, ndim
REAL(dp), DIMENSION(SIZE(p,2)):: psum

CALL amoeba_private

CONTAINS

SUBROUTINE amoeba_private

USE mod_param_initial
USE mod_nrutil

IMPLICIT NONE

INTEGER:: i, ilo, inhi, j 
REAL(dp):: rtol, ysave, ytry, ytmp


ndim = assert_eq(SIZE(p,2),SIZE(p,1)-1,SIZE(y)-1,'amoeba')
iter = 0
psum(:)= REAL(SUM(p(:,:), DIM =1),dp)


!******************************************************************************************
!(1) Iteration loop
!******************************************************************************************
DO

!******************************************************************************************
!(2) Determine which point is the highest (worst), next-highest, and lowest (best)
!******************************************************************************************
  ilo = iminloc(y(:))	!lowest (best) point
  ihi = imaxloc(y(:))	!highest (worst) point 
  ytmp = y(ihi)         !in ytmp the worst value
  y(ihi) = y(ilo)       !in the worst the best
  inhi = imaxloc(y(:))  !second worst (highest)
  y(ihi) = ytmp
  rtol = 2.0_dp*DABS(y(ihi)-y(ilo))/(DABS(y(ihi))+DABS(y(ilo))+tiny)



!******************************************************************************************
!(3) Compute the fractional range from highest to lowest (return if satisfactory)
!******************************************************************************************
!-------------------------------------------------------
!NOTE: if returning, put best point and value in slot 1
!-------------------------------------------------------
	IF (rtol<ftol) THEN 			  
	  CALL swap1(y(1),y(ilo))
	  CALL swap(p(1,:),p(ilo,:))
	  RETURN
	END IF

	IF (iter>=itmax) CALL nrerror('itmax exceeded in amoeba')



!******************************************************************************************
!(4) New iteration
!******************************************************************************************
!--------------------------------------------------------------------------
!(i) extrapolate by a factor -1 through the face of the simplex across 
!    from the worst point, i.e., reflect the simplex from the high point
!--------------------------------------------------------------------------
	ytry = amotry(-1.0_dp)
	iter = iter + 1
    
	OPEN(unit=777, file='track_like.txt',position="append")
	WRITE(777,*)  "Extrapolation by -1", ytry
	WRITE(777,*)  "Iter (extrap -1) = ", iter
	CLOSE(777)
    
	WRITE(3000,*) ""
	WRITE(3000,*) "Iter (extrap -1) = ", iter
	WRITE(3100,*) "Iter (extrap -1) = ", iter

	!OPEN(unit=77, file='track_param.txt', position="append")
	!WRITE(77,*) "Extrapolation by -1", ptry
	!WRITE(77,*) "Iter (extrap -1) = ", iter
	!CLOSE(unit=77)


!--------------------------------------------------------------------------
!(ii) if the result is better than the best point, try an additional 
!     extrapolation by a factor of 2
!--------------------------------------------------------------------------
	IF (ytry <= y(ilo)) THEN 		  		
	  ytry = amotry(2.0_dp)
	  iter = iter + 1
    
	OPEN(unit=777, file='track_like.txt',position="append")
	WRITE(777,*)  "Extrapolation by 2", ytry
	WRITE(777,*)  "Iter (extrap 2) = ", iter
	CLOSE(777)
    
	WRITE(3000,*) "Iter (extrap 2) = ", iter
	WRITE(3100,*) "Iter (extrap 2) = ", iter



!--------------------------------------------------------------------------
!(iii) if the reflected point is worse than the second highest, look for 
!      an intermediate point, i.e., do a one dimensional contraction 
!      (amotry(.) is the likelihood function)
!--------------------------------------------------------------------------
	ELSE IF (ytry >= y(inhi)) THEN  
	  ysave = y(ihi)
	  ytry = amotry(0.5_dp)
	  iter = iter+1


    OPEN(unit=777, file='track_like.txt',position="append")
	WRITE(777,*)  "Extrapolation by 0.5", ytry
	WRITE(777,*)  "Iter (extrap 0.5) = ", iter
	CLOSE(777)
    
	WRITE(3000,*) "Iter (extrap 0.5) = ", iter
	WRITE(3100,*) "Iter (extrap 0.5) = ", iter


!--------------------------------------------------------------------------
!(iv) if this attempt is unsuccessful, contract around the lowest point
!--------------------------------------------------------------------------
   IF(ytry>=ysave) THEN 
	p(:,:)=0.5_dp*(p(:,:)+SPREAD(p(ilo,:),1,SIZE(p,1)))

	DO i = 1, ndim+1
	 IF (i /= ilo) THEN 
	  CALL likelihood_function(p(i,:),y(i),logfct(:))
							
	  OPEN(unit=777,file='track_like.txt',position="append")
	  WRITE(777,*) "Contract likelihood around the lowest point, Row = ", i
	 !WRITE(777,'(F16.4)')   (y(i))
          WRITE(777,*)   (y(i))
	  WRITE(777,*)   "" 
      CLOSE(777)

	  !OPEN(unit=77, file='track_param.txt',position="append")
	  !WRITE(77,*) "Contract likelihood around the lowest point, Row = ",i
      !WRITE(77,'(22F16.4)')   (p(i,j),j = 1, ndim)
      !CLOSE(77)
	  
	 END IF		
    END DO


!--------------------------------------------------------------------------
!(v) keep track of function evaluations
!--------------------------------------------------------------------------
	iter = iter + ndim	
	OPEN(unit=777,file='track_like.txt',position="append")
	WRITE(777,*) "Iter = ", i
	WRITE(777,*)  " " 
    CLOSE(777)
	
	WRITE(3000,*) "Iter = ", i
	WRITE(3100,*) "Iter = ", i			 
	
	psum(:) = REAL(SUM(p(:,:), DIM=1),dp)
    
	END IF
   END IF

 END DO		 

END SUBROUTINE amoeba_private



!******************************************************************************************
!(4) Go back for the test for doneness and the next iteration
!******************************************************************************************
! NOTE: the following function extrapolates by a factor FAC through the face of the simplex
!       across from the high point, tries it, and replaces the high point if the new point 
!       is better

FUNCTION amotry(fac)
   
USE mod_param_initial
USE mod_parameters
USE mod_nrutil

IMPLICIT NONE
!I CHANGED ytry to ytry(1)
REAL(dp), INTENT(IN):: fac
REAL(dp):: amotry
REAL(dp):: fac1, fac2, ytry(1)
REAL(dp), DIMENSION(SIZE(p,2)):: ptry

fac1 = (1.0_dp-fac)/ndim
fac2 = fac1-fac

ptry(:) = psum(:)*fac1 - p(ihi,:)*fac2


!--------------------------------------------------------------------------
!(i) evaluate the function at the trial point
!--------------------------------------------------------------------------
   CALL likelihood_function(ptry,ytry(1),logfct(:))		 
   
   !OPEN(unit=777,file='track_like.txt',position="append")
   !WRITE(777,*) "Extrapolation by a factor of ", fac,":" 
   !WRITE(777,'(F16.4)')    (ytry)
   !WRITE(777,*)  " "
   !CLOSE(777)
   

   !OPEN(unit=77, file='track_param.txt',position="append")
   !WRITE(77,*) "Extrapolation by a factor of ", fac,":"
   !WRITE(77,'(22F16.4)')   (ptry)
   !WRITE(77,*) " "
   !CLOSE(77)


!--------------------------------------------------------------------------
!(ii) if it's better than the highest, then replace the highest
!--------------------------------------------------------------------------
IF (ytry(1)<y(ihi)) THEN	 
  y(ihi) = ytry(1)
  psum(:) = psum(:) - p(ihi,:) + ptry(:)
  p(ihi,:) = ptry(:)
END IF

amotry = ytry(1)

END FUNCTION amotry

END SUBROUTINE amoeba






		